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SUMMARY 


A  new  mixed  finite  element  formulation  is  developed  based  on  the 
Hellinger-Reissner  principle  with  independent  strain.  By  dividing  the  assumed 
strain  into  the  lower  order  part  and  the  higher  order  part,  the  new  formulation 
can  be  made  much  more  efficient  than  the  standard  mixed  formulation.  In  addi¬ 
tion  the  present  new  approach  provides  an  alternative  way  of  introducing  stabi¬ 
lization  matrix  to  suppress  undesirable  kinematic  modes. 
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1.  INTRODUCTION 


Hybrid  and  mixed  finite  element  formulations  based  on  the 
Hellinger-Reissner  principle  or  related  variational  principles  have  been 
available  since  the  early  days  in  the  history  of  finite  element  method  [1].  In 
these  formulations  an  independent  stress  or  strain  is  assumed  within  an  element 
in  addition  to  the  usual  assumed  displacement  expressed  in  terms  of  nodal 
displacement.  The  independent  stress  or  strain  is  eliminated  at  element  level, 
resulting  in  an  element  stiffness  matrix  corresponding  to  nodal  displacement 
vector.  Note  that  in  this  paper  we  are  not  interested  in  the  type  of  mixed  for¬ 
mulation  used  in  reference  9  where  both  nodal  displacements  and  nodal  stress 
variables  remain  in  the  assembled  global  model.  For  an  element  with  a  given 

number  of  nodes  there  is  a  degree  of  flexibility  in  the  choice  of  assumed  stress 

or  strain.  For  example  in  the  case  of  Plan's  original  assumed  stress  hybrid 
model  [2],  the  assumed  stress  is  chosen  to  satisfy  equilibrium  within  an  ele¬ 
ment.  Of  course,  the  property  of  the  stiffness  matrix  depends  very  much  on  the 
assumed  stress  or  strain.  Thus  with  a  proper  choice  of  assumed  stress  or  strain 
it  is  possible  to  develop  a  finite  element  model  which  is  superior  to  the  con¬ 
ventional  finite  element  model  based  purely  on  the  assumed  displacement 
approach.  For  example,  for  thin  plates  and  shells,  a  mixed  formulation  based  on 
the  Hellinger-Reissner  principle  or  the  modified  Hellinger-Reissner  principle 
can  be  used  to  alleviate  the  undesirable  locking  effect  associated  with  the  con¬ 
dition  of  zero  inplane  strain  and  zero  transverse  shear  strain  imposed  on  a 

finite  element  model  f3,4|.  The  role  of  hybrid  and  mixed  formulations  in  con¬ 
junction  with  nearly  incompressible  materials  has  also  been  studied  [5J. 

Moreover  a  mixed  formulation  provides  a  rational  mathematical  basis  for  the 
popular  reduced  and  selective  integration  scheme  f6-8|. 

However,  in  hybrid  and  mixed  formulations  it  is  necessary  to  invert  a 
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matrix  In  order  to  generate  an  element  stiffness  matrix.  Therefore,  In  com¬ 
parison  with  the  assumed  displacement  model  based  on  the  principle  of  virtual 
work,  hybrid  and  mixed  formulations  require  usually  more  time  to  compute  element 
stiffness  matrix  of  the  same  size.  This  fact  has  been  regarded  as  the  major 
drawback  of  the  conventional  hybrid  and  mixed  formulations. 

With  this  in  mind,  we  present  In  this  paper  a  new  mixed  formulation  which 
requires  much  less  computing  time  than  the  conventional  mixed  formulation  for 
the  generation  of  element  stiffness  matrix.  The  new  formulation  is  based  on  the 
Hellinger-Relssner  principle  with  independent  strain.  The  new  formulation  can 
be  applied  to  any  type  of  problems  in  solid  and  structural  mechanics.  However, 
for  simplicity,  two  dimensional  plane  problems  will  be  used  to  demonstrate  the 
effectiveness  of  the  approach.  Initially  nine  node  element  will  be  used  for  the 
purpose  of  illustration.  A  short  discussion  on  four  node  element  will  follow. 

In  order  to  contrast  the  new  formulation  with  the  conventional  formulation,  we 
will  start  off  with  the  conventional  formulation. 


2.  CONVENTIONAL  MIXED  FORMULATION 
2.1  Finite  Element  Formulation 

For  two  dimensional  plane  problem,  the  functional  for  the 
Hellinger-Relssner  principle  can  be  written  as 


r-  .v 


where 


ttd  =  f  (e^  C  e  -  C  e)  dA  -  W 
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(1) 


I  ^xy 


Independent  strain  vector 
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strain  vector  in 
terms  of  displace 
ment  field 


u,  v:  displacement  in  x  and  y  direction  respectively 
A:  area 

W:  external  load  term 

C:  3x3  matrix  of  elastic  constants  integrated  through  thickness 


Note  that  in  eg.  (1)  instead  of  stress,  strain  e  appears  as  independent 
variables  in  addition  to  displacement  field  u  and  v. 


For  finite  element  approximation,  the  displacement  vector  u  is  assumed 
in  terms  of  nodal  displacement  vector  as 


=  NSe 


(2 


where 


N:  shape  function  matrix 


3^:  element  nodal  displacement  vector 


Then  the  strain  vector  c  can  be  written  symbolically  as 


e  = 


(3 


where  B  is  the  matrix  relating  e  to  g.. 


In  addition,  within  each  element  the  independent  strain  e  is  assumed  in  terms  of 
unknown  coefficients  as 


where 


(4) 


P:  strain  shape  function  matrices 

a:  vector  of  unknown  coefficients  in  an  element 


Introducing  eq.  (3)  and  (4)  into  eq.  (1). 

’'r  =  ae  - 


with 


G  =  f  C  B  dA 

~  - 

H  =  f  P^  C  P  dA 

^  ft  ^ 


(5a) 

(5b) 


W  =  2e 


(5c) 


and  Ag  is  the  element  area. 

The  summation  sign  indicates  summation  or  assembly  over  all  elements. 
Taking  =  0  with  respect  a  in  each  element. 


G  <lp  -  H  a  =  0 

A#  *s* 

or  solving  for  a 


a  =  H'^G  qg 


(6) 


(7) 


for  each  element. 


Substituting  eq.  (7)  into  eq.  (5),  is  written  as 

"R  ■  ):(^  E,,  3e  -  al  9e> 

Where  K  =  H’^G 

M  M 

is  the  element  stiffness  matrix. 

Assembling  over  all  elements, 

’'R  ■  '2‘  ~  3  ■  3^9 

where  K:  global  stiffness  matrix 

£:  global  nodal  displacement  vector 
Q:  global  load  vector 

Setting  6iT[^  =  0  with  respect  to  £  leads  to 

<3=9 


(8) 

(9) 

(10) 


which  can  be  solved  for  £,  With  £  thus  £g  known,  the  independent  strain 
vector  e  in  an  element  is  determined  by  substituting  eq.  (7)  to  eq.  (9) 
such  that 

£  =  E!i''S9,  (12) 

and, for  isotropic  materi als, stress  a  is  determined  as  follows: 

a  =  i  C  c  ( 13) 


where  t  is  the  thickness. 


2.2  Nine  Node  Element 

The  present  nine  node  element  is  important  in  that  it  constitutes 
inplane  part  of  the  nine  node  shell  element  described  in  reference  4.  This 
particular  shell  element  has  been  found  to  be  free  of  locking  and  any  com¬ 
patible  or  commutable  kinematic  mode.  Of  course  a  good  nine  node  ele¬ 
ment  is  useful  for  plane  stress  and  plane  strain  analysis  itself.  For  the 
assumed  displacement,  the  element  shown  in  fig.  1  adopts  the  isoparametric 
representation.  As  for  the  coordinate  system,  global  and  local  cartesian 
coordinate  systems  are  used.  The  strains  c,  7  and  displacement  u,v  in  eq. 

(1)  are  defined  with  respect  to  the  global  coordinate  system.  Local  coor¬ 
dinate  system  is  introduced  to  allow  an  assumed  strain  field  with  nonsymmetric 
polynomial  terms.  Local  coordinate  system  is  defined  first  by 
determining  two  unit  vectors  Vj  and  V2  at  c  =  n  =  0  point  such  that 

V  M 

~1  3c  /  3C 


3x  /  3X 

V  =  — -  /  — 

~2  3n '  3n 

where  x  is  the  position  vector  with  components  in  the  global  coordinate 
system.  The  angle  9^  between  these  two  unit  vectors  is  calculated  from  the 
following  equation: 

cose^  =  Vj  •  V2  (15) 
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Then,  if  e^  is  less  than  or  equal  to  90°,  unit  vector  aj  in  the  x  direction 
of  local  coordinate  system  is  defined  as 


35 

H 


Otherwise  a,  is  defined  as 


~1 


3x  /  ax 
3n  f  3n 


(16) 


(17) 


Unit  vector  in  the  y  direction  of  local  coordinate  system  is  orthogonal 
to  a^.  Note  that,  while  Vj  and  V2  determined  at  5  =  n  =  0  point,  a^ 

and  £2  can  be  computed  at  any  point  in  the  element.  Especially  a^  and  ^2 
are  needed  at  integration  points  to  establish  a  local  orthogonal  coordinate 
system.  Now,  for  the  present  nine  node  element,  we  may  assume  the  indepen¬ 
dent  strain  in  local  coordinate  system  as  follows; 

Sx  =  “1  ^  “2^  ^  “S’!  ^  ^  “13^ 

^yy  ®5  “7^  “8^’’  “l4^y 

Sy  =  “9  ^  “10«  ^  ®11^  ^  “12^’’ 

In  eq.  (18),  f^  =  and  f^  =  c^n  for  defined  as  in  eq.  (16).  For  a,^ 

fined  as  in  eq.  (17),  f  and  f  are  chosen  as  f  =  and  f  =  Cn^.  Note 

A  y  X  y 

that,  due  to  the  and  terms,  the  assumed  strains  are  not  symmetric  in  5 

and  n.  However  with  the  use  of  local  coordinate  system  as  defined  here,  the 
resulting  stiffness  matrix  will  not  be  dependent  upon  the  choice  of  global  coor 
dinate  system.  Symbolically  eq.  (18)  can  be  written  as 


e*  =  P*  a 


(19) 
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where 


(19a) 


As  shown  in  eq.  (9),  for  the  generation  of  element  stiffness  matrix,  it  is 
necessary  to  compute  G  and  For  the  present  nine  node  element  the  size 

of  G  and  H  matrices  are  14  x  18  and  14  x  14  respectively  and  G  and  H  matrices 
are  evaluated  using  3x3  point  Gaussian  quadrature.  As  far  as  computung 
time  is  concerned  the  need  to  evaluate  G,  H  and  makes  the  conventional 
mixed  formulation  less  attractive  than  the  assumed  displacement  formulation 
based  on  the  principle  of  virtual  work. 

At  this  point,  it  is  well  to  mention  that  the  and  terms  in  eq. 
(18)  were  chosen  carefully  to  suppress  undesirable  kinematic  modes. 

Without  the  and  terms,  the  assumed  strain  field  is  symmetric  in  c 
and  n.  Moreover,  for  an  element  of  rectangular  shape,  2x2  point  rule  is 
sufficient  enough  for  exact  integration  of  G  and  H.  Then  the  resulting 


stiffness  matrix  will  be  equivalent  to  the  conventional  assumed  displace¬ 
ment  model  with  2  x  2  reduced  integration.  The  equivalence  between  mixed 
formulation  element  and  reduced  or  selective  integration  element  can  be 
proved  by  establishing 

e  =  I  (21) 

at  integration  points.  This  can  be  done  if  the  number  of  integration 
points  is  the  same  as  the  number  of  strain  parameters  or  coefficients  in 
each  component  of  assumed  strain  [3,6,7].  The  element  without  the  s^d 
terms  exhibits  three  kinematic  modes  or  spurious  zero  strain  enery 
modes.  For  a  square  element  with  sides  along  x  =  ±1  and  y  =  ±1  lines, 
these  spurious  modes  are  as  follows  [10,11]; 

(1)  u  =  Cjxd  -  3yd 

V  =  -Cjy(l  -  3xd  (22) 

(2)  u  =  C2(x^  +  y"  -  3xV) 

(3)  V  =  C3(x='  +  y^  -  3xV) 

where  C2  and  C3  are  arbitrary  constants.  Among  these  three  modes,  the 
first  mode  is  incompatible  or  non-commutable  as  it  disappears  for  an 
assembly  of  only  two  elements.  The  second  and  third  modes  are  compatible 
or  commutable  and  may  persist  even  after  assembly  of  elements,  resulting  in 
an  unstable  finite  element  model.  The  and  terms  in  the  assumed 
strain  are  introduced  to  suppress  these  two  compatible  modes.  The  first 
mode  is  left  in  the  element  since  it  is  harmless.  It  is  to  be  noted  that 
the  same  assumed  strain  was  used  in  the  shell  element  described  in 
reference  4  where  it  was  found  to  be  very  effective  in  alleviating  locking 
effect  associated  with  the  condition  of  zero  inplane  strain. 
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3.  NEW  FORMULATION 


The  proposed  new  formulation  is  also  based  on  the  Hellinger-Reissner 
principle  with  the  functional  given  in  eq.  (1).  However  in  the  new  for¬ 
mulation,  the  independent  strain  is  written  as 

£  =  £l  +  £h  (23) 

where  is  the  independent  strain  vector  with  lower  order  polynomial  terms 
in  5  and  n.  On  the  other  hand  e^,  the  higher  order  strain  vector,  contains 
higher  order  terms  in  c  and  n.  More  discussion  on  and  will  be  given 
later.  Inserting  eq.  (23)  into  eq.  (1),  the  functional  Tr|^  becomes 

’'r  =  /(£l  £  I  - 

'•■  /  £h  £  £  /  £H  ~ 


-/7£5££HdA.W  (24) 

To  help  illustrate  the  new  formulation,  we  will  consider  again  nine  node 
element.  For  nine  node  element  of  rectangular  shape,  the  highest  order 
terms  in  the  7  matrix  are  quadratic  in  either  c  or  n  direction.  Thus  if  the 
matrix  is  bilinear,  then  the  first  integral  can  be  integrated  exactly 
using  2x2  point  Gaussian  quadrature.  For  e^,  we  start  by  assuming  ejj, 
the  higher  order  strain  vector  in  the  local  coordinate  system,  as 


(Sx^H  "  “l^x 
(eyy)H  =  «2^y 


(25) 
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Or  in  matrix  form. 


*  _* 

Em  =  P  a 

#vn  ^  ^ 


where 


(26) 


(26a) 


ot2J  (26b) 


ic 

Note  that  contains  the  higher  order  and  terms  of  the  conventional 
formulation  given  in  the  previous  section.  Using  the  strain  transformation 
matrix  T,  the  vector  in  the  global  coordinate  system  is  expressed  as 


^  I  Sh 


=  T  /a 


(27) 


SVbVViV 
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.r-WiV 


where  "P  =  T  P  (27a) 

Of  course,  for  rectangular  element,  7  is  equal  to  7.  However,  for  an  ele¬ 
ment  of  arbitrary  shape,  strain  transformation  in  eq.  (27)  is  used.  Again, 
for  rectangular  element,  the  third  integrals  in  eq.  (24)  require  2x2 
point  rule.  On  the  other  hand  the  second  and  fourth  integrals  require  3  x 
3  point  rule.  For  an  element  with  arbitrary  shape  the  argument  regarding 
the  number  of  integration  points  for  exact  integration  does  not  hold. 
However  even  in  this  case  the  same  integration  rules  will  be  adopted. 
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In  addition,  for  the  first  and  third  integrals,  we  may  write  £|^  as 


Jl  =  J^n^(5,n)e, 


(28) 


where  are  bilinear  shape  function  such  that  N.  =  1  at  the  point  i  of  the 
2x2  point  integration  and  zero  otherwise,  and  1.  is  the  strain  7  evaluated 
at  integration  point  i.  In  another  word,  we  can  set 


El  -  e 

at  2  X  2  integration  points. 

Introducing  eq.  (29)  to  eq.  (24), 


(29) 


TTn  = 


/  -ley  C  Cl  dA  +  f  Cu  C  cdA 

2-vL  ~  ~L  ~H  ~  ~ 


eu  C  Ci  dA 

•srn  ^ 


(30) 


In  eq.  (30),  letters  L  and  H  under  the  integral  signs  indicate  lower  order 
integration  (2x2  point  rule)  and  higher  order  integration  (3x3  point 
rule)  respectively.  Introducing  eqs.  (3)  and  (26)  into  eq.  (30)  and  noting 
eq.  (29),  can  be  written  as 


(31) 


where 


K,  =  r  B^C  B  I  J  I  dc  dn 

^  ^  /V  /V  I  ^  I 


6  =  Gii  ”  G| 
““  ~n  "'L 


(31a) 

(31b) 
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r  p^c  B  I  j 

^  ^  M  I 

n 

/  P^C  B  I  J 

^  M  ^  ^  I 

/  "P^C  T^  I  J 

^  ~  ~  ~  ~ 


dc  dn 


dc  dn 


dc  dn 


(31c) 

Old) 

(39e) 


and  I  J  I  is  the  determinant  of  Jacobian  matrix  J.  Note  that,  although  the 
same  symbol  B  appears  in  K^,  and  Sl*  &  evaluated  at  3  x  3 

Gaussian  integration  points  while  B  in  and  6^^  is  evaluated  at  2  x  2 

Gaussian  integration  points.  However  to  save  computing  time,  the  B  matrix 
at  2  X  2  integration  points  can  be  interpolated  from  the  B  matrix  evaluated 
at  3  X  3  integration  points.  That  is,  we  evaluate  B  at  2  x  2  points  from  the 
following  expression: 


B  * 


?(i(5.Ti)Bi 


(32) 


where  B^  is  the  B  matrix  at  the  integration  point  i  and  is  the  shape 
function  such  that  N^-  =  1  at  point  i  of  the  3x3  integration  points  and 
zero  at  other  points.  In  addition  the  determinant  |  J|  at  2  x  2  integration 
points  is  also  interpolated  from  |  J  |  at  3  x  3  integration  points. 

Taking  5Tr|^  with  respect  to  a  of  each  element 

H  2  -  &e  =  0  (33) 


or  solving  for  a 

2  ■  r'  I  a.  (S'!) 

for  each  element. 
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Substituting  eq.  (34)  into  eq.  (31), 


-R  =  Ki  K,  Se  -  al  9,)  (35) 

where 

!le  ■  Kl  I  <36) 

is  the  element  stiffness  matrix.  Assembling  over  all  elements 

•r  ■  3  -  3^3  (33) 

where  K  is  the  global  stiffness  matrix  corresponding  to  the  global  displa¬ 
cement  vector  3  and  Q  is  the  global  load  vector.  Taking  6^^  =  0  with 
respect  to  ^  results  in 

K  a  •  3  (38) 

which  can  be  solved  for  £.  With  £  and  thus  £g  known  strain  e  is  determined 
from  eqs.  (23),  (27),  (28)  and  (34)  as  follows: 

e  =  Cl  +  Cu  =  El  £® 

~L  ~L  ~~ 


Nj  B,  a.*£K'‘S3e 


Then  stress  a  is  determined  by  eq.  (13). 

It  should  be  pointed  out  that  the  element  stiffness  matrix  in  eq.  (36) 

has  two  components,  K|^  and  element  of  rectangular  shape  the 

K.  matrix  is  the  same  as  the  stiffness  matrix  of  the  conventional  assumed 
~L 

displacement  model  with  2x2  reduced  integration.  Again  the  equivalence  is 
established  through  eq.  (29).  Then  the  matrix  has  the  same  three  kinematic 


modes  given  in  eq.  (22).  With  the  addition  of  ^matrix,  the  compatible 

kinematic  modes  are  suppressed,  leaving  only  the  incompatible  mode.  Therefore, 
as  far  as  kinematic  modes  are  concerned,  both  the  conventional  formulation  and 
the  new  formulation  result  in  the  same  incompatible  mode.  However  element 
stiffness  matrix  from  the  new  formulation  is  not  exactly  the  same  as  that  from 
the  conventional  formulation.  The  size  of  ^  and  F  matrices  in  the  new  for- 
mulation  are  2  x  18  and  2x2  respectively.  Therefore,  computation  of  the 
6^  H"^  G  matrix  in  the  new  formulation  can  be  carried  out  without  much  effort. 

A#  ^ 

Recall  that  in  the  conventional  formulation,  the  sizes  of  G  and  H  matrices  are 
14  X  18  and  14  x  14  respectively.  Note  also  that  the  present  element  passes  the 
patch  test. 

4.  COMPARATIVE  NUMERICAL  TEST 

(a)  Comparison  of  Computing  Time 

In  order  to  evaluate  computing  efficiency  of  the  new  formulation,  a 
test  was  run  in  which  stiffness  matrix  of  single  nine  node  element  was  com¬ 
puted  40  times  consecutively.  Table  1  shows  relative  computing  time  for 
different  element  types.  Clearly  the  new  mixed  formulation  element  (NM) 
requires  much  less  computing  time  than  the  conventional  mixed  formulation  ele¬ 
ment  (CM).  Surprisingly  the  new  mixed  formulation  takes  less  computing 
time  than  the  conventional  assumed  displacement  element  with  3x3  point 
rule  (DISP3).  Of  course  the  assumed  displacement  element  with  2x2  point 
rule  (DISP2)  takes  the  least  time.  However,  as  mentioned  before,  this  element 
has  compatible  kinematic  modes  and  thus  cannot  be  used  in  general  stress 
analysis. 

(b)  A  Panel  under  a  Horizontal  Point  Load 

Figure  2  shows  a  rectangular  panel  subjected  to  a  horizontal  point  load  P. 
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The  panel  Is  modeled  by  2  x  12  mesh  as  shown  In  the  figure.  Two  boundary  con¬ 
ditions  are  considered.  In  case  1.  the  left  end  Is  completely  fixed.  In  case 
2,  the  vertical  displacements  of  the  first,  second,  fourth  and  fifth  nodes 
along  the  left  end  are  left  unconstrained.  The  first  case  has  been  used  In 
references  13  and  14  to  demonstrate  the  undesirable  effect  of  spurious  kinematic 
modes  on  numerical  solution.  The  pertinent  data  are  as  follows; 

length  L  «  12m 
depth  b  «  2m 
Poissons  ratio  u  =  0.2 

*  W 

where  E  Is  the  Young's  modulus  and  A  Is  the  cross-sectional  area.  Table  2  lists 
the  nondimenslonal  horizontal  displacement  calculated  at  the  load  point  for  the 
NM,  CM.  DISP2  and  D1SP3  elements.  Numerical  solutions  are  nondimenslonal 1  zed  by 
dividing  the  computed  values  by  the  tip  displacement  PL/AE  of  a  panel  under  total 
load  P  distributed  uniformly  over  the  cross-sectional  area.  The  solution  for  the 
DISP2  element  Is  very  large  compared  with  those  obtained  by  the  other  three  ele¬ 
ment  types,  especially  In  case  2.  This  Is  due  to  the  spurious  compatible  kine¬ 
matic  modes  In  the  01SP2  element  triggered  by  the  point  load.  The  NM,  CM  ele¬ 
ments  are  free  of  compatible  kinematic  modes  and  the  0ISP3  element  has  no 
kinematic  mode.  Therefore  they  provide  stable  solutions. 

(c)  A  Cantilever  Beam 

A  cantilever  beam  subjected  to  a  tip  load  P  1s  used  to  evaluate  the 
performance  of  the  present  new  element  as  compared  to  the  conventional  mixed 
formulation  element  and  the  assumed  displacement  model  element.  Note  again  that 
the  same  problem  has  been  used  as  a  numerical  example  In  reference  13.  As 
Illustrated  In  fig.  3,  the  cantilever  beam  Is  modelled  by  three  different 
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finite  element  meshes.  Two  different  length  to  depth  ratios  of  10  and  20 
are  considered.  Table  3  lists  computed  vertical  displacement  of  point  A  in 
fig,  3  and  stress  evaluated  at  point  B.  The  point  B  is  one  of  the 
Gaussian  points  in  the  2x2  integration  rule.  Numerical  solutions  were 
normalized  by  the  following  solutions  obtained  from  the  Bernoulli -Euler 
beam  theory: 

v^  =  plVsei 
(a^x^B  “  ” 

where 

I  =  sectional  area  moment  of  inertia 
M  s  f'ending  moment 
yg  *  y  coordinate  of  point  B 

Numerical  results  in  Table  3  indicates  that  NM,  CM  and  DISP2  elements  perform 
much  better  than  the  0ISP3  element  especially  for  meshes  with  non-rectangular 
elements.  For  this  particular  example,  the  NM  element  seems  to  be  better  than 
the  CM  element.  It  is  interesting  to  note  that,  for  the  present  problem  under 
vertical  tip  load,  the  spurious  kinematic  modes  of  the  DISP2  elements  remain 
untriggered  as  Indicated  by  the  stable  solution.  This  is  in  sharp  contrast  to 
the  previous  example  under  horzontal  point  load. 

5.  FOUR  NODE  ELEMENT 

For  four  node  element,  a  finite  element  model  based  on  the  conventional 
mixed  formulation  may  be  developed  by  assuming  independent  strain  vector 
e*  defined  with  respect  to  the  local  coordinate  system  as  follows: 

Sx  '  “l  *  “49x 
‘Jy  ■  “Z  *  “59y 
®xy  «  «3 


(42) 


In  eq.  (42),  =  n  and  gy  =  C  for  defined  as  in  eq.  (16).  For  a^  defined  as 

in  eq.  (17),  g  and  g„  are  chosen  as  g^  =  C  and  g„  =  ti.  Again  strain  e  in  the 
global  coordinate  system  is  obtained  from  e*  through  strain  transformation  as 
shown  in  eq.  (20).  As  for  numerical  integration,  the  G  and  H  matrices  in  the 
conventional  formulation  are  evaluated  by  2  x  2  point  rule.  For  a  finite  ele¬ 
ment  model  based  on  the  new  formulation,  the  higher  order  strain  in  the  local 
coordinate  system  is  assumed  as 


“  “l^x 
(®yy)H  “  ®2^y 


Again  in  the  global  coordinate  system  is  determined  from  through 

strain  transformation.  As  for  numerical  integration,  the  and 

G|_  matrices  are  evaluated  by  one  point  integration  rule  whereas  2x2  point 

rule  is  used  for  integration  of  G^  and  ¥  in  the  new  formulation.  Of  course 

in  computing  Kj^  and  B  matrices  and  |  J  |  at  the  integration  point  can  be 

interpolated  from  B  matrix  and  |  J  |  evaluated  at  2  x  2  initegration  points 

following  the  similar  procedure  used  for  the  nine  node  element.  For  an  element  of 

rectangular  geometry,  the  matrix  is  the  same  as  the  stiffness  matrix  of  four 

node  assumed  displacement  element  based  on  the  principle  of  virtual  work  with 

one  point  integration.  The  matrix  has  two  compatible  spurious  kinematic 

modes.  However  addition  of  the  G^H'^g  matrix  as  shown  in  eq.  (36)  suppresses 

these  kinematic  modes  and  thus  element  stiffness  matrix  is  stable  and  has 

~e 

correct  rank.  For  an  element  of  rectangular  shape,  the  new  formulation  element  is 
equivalent  to  the  element  based  on  the  conventional  formulation.  Furthermore, 
it  is  found  that  for  an  element  with  square  geometry,  the  present  element  is 
equal  to  the  Belytschko's  four  node  element  with  a  properly  chosen  stabilization 
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r»i^j  «r*  •*-.  T'.- 1 


matrix  [12j.  The  element  stiffness  matrix  of  Belytschko's  four  node  element 
described  in  reference  12  may  be  expressed  as 

^  ^  ♦  Ks  («) 


where 


!ls 


Et 

1-v' 


r  T 


'2X2X2 


(44) 


is  the  stabilization  matrix  introduced  to  suppress  kinematic  modes  in  Kj^, 

The  constants  Cj  and  C2  are  control  parameters  and  the  expression 
for  and  ^2  given  in  reference  12,  If  we  set  Cj  *  “  1/12,  then 

for  a  square  element,  the  resulting  stiffness  matrix  is  exactly  the  same  as 
that  of  the  present  four  node  element. 

6.  DISCUSSION  AND  CONCLUSION 

Numerical  test  with  nine  node  element  indicates  that  the  proposed  new 
formulation  needs  less  than  half  of  the  computing  time  required  for  the 
conventional  mixed  formulation  to  generate  element  stiffness  matrix.  Also  for 
nine  node  element  the  computing  time  for  new  element  is  slightly  less  than 
that  required  for  the  conventional  assumed  displacement  model  with  3  x  3 
point  integration.  The  nine  node  element  based  on  the  new  formulation  is 
not  exactly  the  same  as  that  based  on  the  conventional  mixed  formulation. 
However  they  are  very  close  to  each  other.  For  four  node  element,  the 
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p: 
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conventional  and  new  formulation  result  in  the  same  stiffness  matrix  for  rec¬ 
tangular  element  geometry.  Also  equivalence  to  Belytschko's  four  node  ele¬ 
ment  with  properly  adjusted  stabilization  control  parameters  has  been 
observed.  It  is  to  be  noted  that  the  ff'^G  matrix  in  eq.  (36)  associated  with 
higher  order  assumed  independent  strain  plays  the  role  of  stabilization  matrix. 
As  such  the  present  new  approach  can  be  viewed  as  an  alternative  way  of  intro¬ 
ducing  stabilization  matrix  into  the  finite  element  formulation.  The  pre¬ 
sent  formulation  can  be  easily  extended  to  two  and  three  dimensional 
problems,  as  well  as  thin  plate  and  shell  problems.  In  fact  a  new  approach 
applied  to  shell  element  formulation  will  be  the  subject  of  a  forthcoming  paper. 
Also  extension  to  nonlinear  problems  seems  to  be  straightforward. 
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Table  2.  Nondimensional  Horizontal  Oispl 
a  Cantilever  Panel. 


Element  Type 

NM 

CM 

Case  1 

1.3087 

1.2456 

Case  2 


1.3096 


1.2466 


at  the  Load  Point  of 


DISP3 

DISP2 

1.2116 

147.62 

I5 


1.2126 


7.4928  X  10 


[•TlTr 


\L 


Mesh  Type 

1 

2 

3 

L/b 

Type 

% 

% 

"a 

% 

% 

NM 

.9950 

1.0223 

1.0141 

1.0848 

.9875 

.9805 

CM 

.9900 

1.0000 

.9748 

.9140 

.9603 

.9228 

10 

DISP3 

.9541 

1.1407 

.7913 

.6957 

.7370 

.7745 

DISP2 

1.0058 

1.0000 

1.1085 

1.1253 

.9549 

.9584 

NM 

.9902 

1.0223 

.9905 

1.0045 

.9833 

.9933 

CM 

.9850 

1.0000 

.9672 

.9003 

.9556 

.9228 

20 

DISP3 

.9362 

1.1931 

.7584 

.6833 

.4406 

.4871 

DISP2 

1.0014 

1.0000 

1.1036 

1.1253 

.9506 

.9584 

